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Grand canonical molecular dynamics (GCMD) is applied to the nucleation process in a metastable 
phase near the spinodal, where nucleation occurs almost instantaneously and is limited to a very 
short time interval. With a variant of Maxwell's demon, proposed by McDonald [Am. J. Phys. 31 
(1963): 31], all nuclei exceeding a specified size are removed. In such a steady-state simulation, 
the nucleation process is sampled over an arbitrary timespan and all properties of the metastable 
state, including the nucleation rate, can be obtained with an increased precision. As an example, 
a series of GCMD simulations with McDonald's demon is carried out for homogeneous vapor to 
liquid nucleation of the truncated-shifted Lennard-Jones (tsLJ) fluid, covering the entire relevant 
temperature range. The results are in agreement with direct non-equilibrium MD simulation in the 
canonical ensemble. It is confirmed for supersaturated vapors of the tsLJ fluid that the classical 
nucleation theory underpredicts the nucleation rate by two orders of magnitude. 

I. INTRODUCTION 

The key properties of nucleation processes are the height Af2* of the free energy barrier that must be overcome to 
form stable embryos of the emerging phase and the nucleation rate J that indicates how many nuclei appear in a given 
volume per time. The most widespread approach for calculating these quantities is the classical nucleation theory 
(Fede r et all Il966). which has significant s hortcomings, e.g., it can overestimate Ai?* significantly for homogeneous 



vapor to liquid nucleation ( Talanquerl . 120071 ). A more accurate theory of homogeneous nucleation, which is sought after, 



would also increase the reliability for more complex applications such as heterogeneous and ion-induced nucleation in 
the earth's atmosphere. 

An important problem of the classical nucleat ion theory (CNT) is that the underlying basic assumptions do not 
apply to n anoscopic nuclei llVrabec et a 1 2009^1 . Although it is possible to measure the critical size by neutron 



scattering (jPebenedettil . 1200a Pan et al. , l2006t ). the thermophysical properties of such nuclei are mostly very hard 



to i nvestigate experimentally. However, they a r e wel l accessible by calculations based on density func tional the- 
ory ijOxtobv and Evansl. Il988l: IZeng and Oxtobvl. Il99ll: iBvkov and S hchek ml Il999l: lUline and Cortil . l2008h as well as 
molecular simu lation ll Vrabec et alll2006t|Hoivst and Litniewskil.l2008tlSchrader et alll2009f). For ins tance, vaporiza- 
tion processes ijHolyst and Litniewskil . l2008l ) and equilibria ( Vrabec et alll2Q0^ Schrader et all l2009l l of single liquid 



droplets can be simulated to obtain the surface tension as well as heat and mass transfer properties of strongly curved 
interfaces. Similarly, very fast nucleation processes that occur in th e immediate vicinity of th e spino dal are experi- 
mentally i naccessible, whereas they can be studied by Monte C arlo ( Neimark and Vishnvakovl . l2005l l and molecular 



dynamics (jYasuoka and Matsumotol . Il998l : iHorsch et all l2008h simulation of systems with a large number of parti- 
cles. Lower nucleation rates are accessible by tr ansition path sampling based methods such as forward flux sampling 
( Ghiringhelli et all . l2008t Ivan Meel et a.1.1 . l2008h . Hence, molecular simulation is crucial for the further development 



of nucleation theory. 

Such molecular dynamics (MD) simulations, dealing with single nuclei in equilibrium as well as with homogeneous 
nucleation processes in supersaturated vapors, led to the formulation of a surfac e propery corrected (SPC) modifi- 
cation of CNT for vapor to liquid nucleation of unpolar fluids, cf. previous work ( Horsch et all 12008?) for a detailed 



presentation and justification. Both CNT and the SPC modification apply the expression 



dQ\ (dT\ 

dv ) p,y,T \ dv ) \,V,T 



accounting for the positive contribution of the surface tension 7 acting on the surface area T as well as a negative 
bulk contribution, where fi and ^ a {T) are the chemical potential of the supersaturated and the saturated vapor at 
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the temperature T, respectively, to the free energy of formation 



An v =i (2) 



for a nucleus containing v particles. The maximal free energy of formation A,!?*, corresponding to the critical size 
v* , is the decisive quantity for the nucleation rate, given by the Arrhenius equation as 

J^exp^j, (3) 

where fee is the Boltzmann constant. In addition to the usual collision term from kinetic gas theory, the pr e-exponentia l 



coeffi cient A includes the Zel'dovic (SejibflOBira) factor and a correction for thermal non-accomodation IjFeder et all 
1966). Although A is not constant, it depends to a much lower extent on supersaturation than the exponential term. 
CNT applies the capillarity approximation, which in the present context means that 7 is assumed to be the same 
as the surface tension of the planar vapor- liquid interface 700 . The surface area is determined from the assumption 
that all nuclei are spherical. The SPC modification replaces the capillarity approximation with the Tolman equation 
(|Tolmanl . [l949l l. 



+ «> 

wherein 1Z is the radius of the nucleus and S is the Tolman length, a characteristic interface thickness, while the 
surface area is increased by a steric factor s. In particular, the temperature-dependent correlations 

with respect to the critical temperature T c , as well as 

_ 0.85 (1 - r/Tc)' 1 + (v/75) 1 / 3 

l + (i//75)V3 ' [b) 

can be used for unpolar fluids 1 Horsch et all 120081 ) . A different appro ach is given by the Hale scaling law (HSL) . In 



can be used tor unpolar nuids l|Horscn et all I^UUal l . A ditterent appro acn is given by tne Male 
agreement with experimental data on nucleation of water and toluene ()Hald . Tl986l ). it predicts 



\27(fi - \i a 

where the proportionality constant only depends on properties of the critical point. 

The present work has the objective of refining the methodology used for direct MD simulation of nucleation pro- 
cesses. According to the method of Yasuoka and Matsumoto (YM), a supersaturated vapor is simulated in the 
canonical ensemble and the nucleation rate is obtained from the number of n uclei formed over time, using a linear 
fit where only nuclei that exceed a sufficiently large threshold size are counted ( Yasuoka and Matsumotol . Il998h . Nu- 



cleation occurs after the metastable state is equilibrated and before nucleus growth becomes dominant. However, 
the timespan corresponding to nucleation is very short for the high nucleation rates that are accessible to direct MD 
simulation, which restricts the statistical basis and the precision of the results. Near the spinodal, the regimes of 
equilibration, nucleation, and growth even start to overlap and the YM method becomes unreliable. 

Wedekind et al. recently developed a more rigorous me thod which is based on mean first passage times (MFPT) 
obtained by averaging over hundreds of simulation runs (|Wedekind et all [2007). But as Chkonia et al. point out, 



'the computational costs of making the necessary repetitions to evaluate the MFPT can be very high,' whereas YM 
requires many clusters f orming and it therefo re becomes more sensible to deviations coming from vapor depletion or 
coalescence of clusters' ( Chkonia et al l l2009l l. 



A new direct equilibrium MD simulation method is introduced in the present work. The underlying concept is to 
simulate the non-equilibrium as a stationary process in the grand canonical ensemble. Thereby, it is possible to sample 
exclusively nucleation as opposed to nucleus growth and coalescence. While the precision of the results is increased by 
maintaining the steady state over an arbitrarily long time interval, the advantages of the YM non-equilibrium method 
are also retained. In particular, only one MD simulation run is required and the nucleation rate is obtained from the 
number of larg e nuclei formed over time. This is achieved by combining grand canonical molecular d ynamics (GCMD) , 
introduced by I CielinskJ . Il985h . and an 'intelligent being' that continuously removes all large nuclei l|McDonaldl .ll962): 



McDonald's demon. 



3 



II. SIMULATION METHOD 

In a closed system, nucleation is an instationary process because the metastable phase is depleted by the emerging 
nuclei. The idea behind the present approach is to simulate the production of nuclei up to a given size for a specified 
metastable state. Nuclei above the given size are extracted, and particles are inserted as monomers into the system 
to replenish the metastable phase. 

GCMD regulates the chemical potential and samples the grand canonical ensemble: alternating with standard 
MD steps, particle s are deleted from and inserted into the system probabilistically with the usual grand canonical 
acceptance criteria fCielin skl , ll985t[Lupkowski and van Swoill99ll ]. For a test deletion, a random particle is removed. 



For a test insertion, the coordinates of an additional particle are chosen at random. The potential energy difference 
5V is determined for each of the test operations and compared with the residual chemical potential. The acceptance 
probability is defined the same way as for the Metropolis algorithm, i.e., it is 



P = min pA exp 



-/i - SV 



k B T 



,1 , (8) 



in case of deletions and similar for insertions (Al len and Tildeslevl . ll987l ). In this expression, p is the density and A 



is the thermal wavelength. The number of test deletions and insertions per simulation time step was chosen in this 
work between 10 -6 and 10 -3 times the number of particles. 

Whenev er a nucleus exceeds th e specified threshold size (9, McDonald's demon (McDonald], [L962) - called Szilard's 
demon by i Schmelzer et al.l . ll997t ) - removes it from the system and replaces it by a representative configuration of the 



metastable phase. If a dense phase is simulated, this can be achieved by, e.g., inserting an equilibrated homogeneous 
configuration in the center of the free volume, followed by preferential test insertions and deletions in the affected 
region. In a supe rsaturated vapor, however, the density is usually so low that it is sufficient to leave a vacuum behind 
as suggested by (|McDonaldl . ll962l l. 

Establishing a steady state by continuously removing the largest nuclei is the purpose and the main advantage of 
McDonald's demon. Consequently, the further behavior of these nuclei cannot be tracked. It is assumed that most 
of the nuclei that are extracted would have continued to grow and that the demon intervention rate Je is therefore 
similar to the actual nucleation rate J. The deviation between these rates can also be quantified by regarding the size 
evolution of a single nucleus in terms of a discrete one-dim ensional random wal k over the order parameter v. At each 
size transition, v is either decreased or increased by one tH orsch et all I2009T ). The short-term growth probability, 



corresponding to a size increase in the next step, is then given by 



where Z v ±\ is the grand canonical partition function under the condition that the nucleus contains v ± 1 particles. 
Neglecting all discrete size effects, the long-term growth probability q u , which corresponds to the ca ses where the 



nucle us never evaporates completely and thus eventually reaches arbitrarily large sizes, has the property ijHorsch et all 
l2009h 

dZ _ -d(dq u /dv) 



Zdv 2 (dq v / dv) dv 

Using adequate boundary conditions, the long-term growth probability of the nuclei that are removed by the demon 
can be determined as 

Ji Z- 2 dv 

le = r oo „_ 2 , • (11) 
J 1 Z v dv 

The intervention rate is therefore related to the nucleation rate by 

f /2AQ U \ , f°° f2AO v \ , , . 

Je k ^{-^T) dv=J L cxp u^r- (i2) 



In particular, for a threshold size sufficiently above v* the approximation J w Je is valid ( Horsch et all [20091. 

The truncat ed-shifted Lennard- Jones (tsLJ) fluid accurately describes the fluid phase coexistence of noble gases 
and methane i Vrabec et all l2006t ). avoiding long-range corrections which are tedious for inhomogeneous systems. 



Homogeneous vapor to liquid nucleation of the tsLJ fluid was studied here by GCMD simulation with McDonald's 
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demon at temperatures of 0.65 to 0.95 in units of e/fce (where e is the energy parameter of the Lennard- Jones 
potential). Note that the triple point t emperature of the tsLJ fluid is Ta = 0.65 w hile T c is 1.078 so that the entire 
relevant temperature range is covered ( Vrabec et al.l . 120061 : Ivan Meel at, al.l .l2008h. The ((Stmingei [l963) criterion 
was used to discern the emerging liquid from the surrounding supersaturated vapor and nuclei were determined as 
biconnected components. 



III. SIMULATION RESULTS 

Figure Q] shows the aggregated number of demon interventions in one of the present GCMD simulations and, for 
comparison, the number of nuclei in a MD simulation of the canonical ensemble under similar conditions. The constant 
value of the supersaturation 

S = CXP { k B T J' (13) 

in the GCMD simulation agreed approximately with the time-dependent S in the NVT simulation about t = 400 
after simulation onset in units of a(m/e) 1 / 2 , wherein a is the size parameter of the Lennard- Jones potential and m is 
the mass of a particle. 

During the NVT run, however, S decreased from about 3 to 1.5. The observed rate of formation was significantly 
lower for larger nuclei, which is partly due to the the depletion of the vapor over simulation time. Depletion causes less 
monomers to interact with a nucleus surface when large nuclei are formed because by that time, a substantial amount 
of particles already belong to the liquid. Moreover, a small nucleus will eventually decay with a higher probability, 
given by 1 — q, instead of growing to arbitrarily large sizes, cf. Eq. ifTTj) . Therefore, large nuclei are necessarily formed 
at a lower rate. 

In Fig. it can be seen how the decreasing supersaturation in the canonical ensemble MD simulation affects the 
nucleus size distribution. Around t = 400, the distribution of small nuclei present per volume was similar in both 
simulation approaches. Near and above the critical size, i.e., 27 particles according to CNT, cf. Tab.[H deviations arise 
because of the different boundary conditions. Comparing the distribution for the grand canonical steady state with the 
corresponding theoretical predi ction shows that CNT underestimates the number of nuclei present in the metastable 
state, confirming the result of (jTalanquerl . l2007l l that CNT exaggerates the fre e energy of nucleus formation. 

CNT is also known to underestimate the nucleation rate of unpolar fluids i Horsch et all 120081 ). The determined 



demon intervention rates confirm this conclusion, cf. Tab. [J and as shown in Fig. [3j the HSL is significantly more 
accurate than CNT for low temperatures. For T = 0.85, HSL and CNT lead to similar predictions, deviating from 
simulation results by two orders of magnitude. At T = 0.95, a nucleation rate of In J = -16.08 was obtained for 
S = 1.146 (using > 3 i/g PC ) where CNT predicts lnJcNT = -19.99, cf. Tab. U as opposed to lnJksL = -24.27. 
Thus, HSL breaks down at high temperatures for the tsLJ fluid. Present results generally agree with nucleation rates 
obtained by NVT simulation at temperatures between 0.65 and 0.95, as can be seen by comparison with the SPC 
modification that was correlated to data from canonical ensemble MD simulation ( Horsch et all 12008). 



Figure |4] shows how the choice of affects the nucleus temperature. The largest nuclei allowed to remain in the 
system have a highly elevated temperature and the amount of nucleus overheating can be explained by considering 
the boundary condition that McDonald's demon imposes on size fluctuations. Only nuclei that do not fluctuate to 
sizes above remain in the system for a significant time. Almost all nuclei with v k. approach the point where 
overheating due to the enthalpy of vaporization released during condensation countervails the supercooling of the 
vapor. Note that this effect is much stron ger than the overheating AT* of the critical nucleus according to CNT due 
to nucleation kinetics ( Feder et all fl966) 



AT* = 2 -^, (14) 
Ah v v ; 

where fz is the ZePdoviC factor and Ah v is the enthalpy of vaporization, evaluating to AT^ NT = 0.00608 in the 
present case. 

With a threshold far below the critical size, the intervention rate of McDonald's demon is several orders of magnitude 
higher than the steady-state nucleation rate, cf. Tab. Inland Fig.[H In agreement with Eq. lfl2"]) . Jq reaches a plateau 
for > i/g pc . In particular, the approximation J ss Jq is valid for all values shown in Tab. U and Fig. [31 As Tab. 
Ull also shows, the density and the pressure of the supersaturated vapor have very good convergence properties with 
respect to the intervention threshold size and can already be accurately obtained at a high accuracy for values near 
the critical size. 
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IV. CONCLUSION 

GCMD with McDonald's demon was established as a method for steady-state simulation of nucleation processes. 
The main purpose of the new method consists in directly simulating a metastable state that undergoes a phase 
transition at a high rate without being limited to sampling only the short timespan until nucleation occurs. 

By implication, growth or decay processes of very large nuclei are not covered. These hava to be considered using 
the cutoff correction given by Eq. (fT2| unless the intervention threshold size is significantly larger than v*. Due 
to an intervention scheme based on the single order parameter u, other relevant order parameters such as shape or 
temperature of the nuclei can experience a perturbation for a nucleus size similar to 0. It was shown for the nucleus 
temperature that this only concerns the largest nuclei in the system and that the range of nucleus sizes unaffected by 
intervention based overheating can be extended arbitrarily if a sufficiently high value of is chosen. 

The intervention rate necessarily approaches the nucleation rate for increasing values of the intervention threshold 
size. The dependence of Jq on is already accurately described for > z/*/2 by modeling the nucleus size evolution 
as a one-dimensional random walk without taking any other order parameter into account. 

For vapor to liquid nucleation of the tsLJ fluid, a series of simulations was conducted over a wide range of tem- 
peratures. Good agreement with canonical ensemble MD simulation results was reached. It was confirmed that 
CNT overstates the free energy of nucleus formation and underpredicts the nucleation rate. HSL accurately describes 
nucleation near the triple point temperature; at high temperatures, however, significant deviations are present. 
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Table I 

Average number of particles and intervention rate of McDonald's demon during GCMD simulation as well as the 
nucleation rate approximated by J « 9©(CNT)Je, cf. Eqs. ifTTj) and (jT2j) , in dependence of simulation conditions, 
i.e., temperature (in units of e/fe), supersaturation, intervention threshold size (in particles), and system volume (in 
units of cr 3 ), compared to theoretical predictions for the nucleation rate according to CNT and the SPC modification; 
all logarithms are given with respect to the reduced rates, normalized by (m/e) 1 / 2 ^ 4 . Note that the intervention 
threshold size is sufficiently larger than the critical size in all cases. 
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Table II 

Pressure supersaturation p/p a {T) and density supersaturation p/p a (T) as well as the intervention rate of McDonald's 
demon in dependence of simulation conditions along with the long-term growth probability qe of a nucleus containing 
particles, cf. Eq. ifTTj). according to CNT. 
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Figure Q] Top: number per unit volume p n of nuclei containing v > 25 (• • -), 50 (— ), and 150 ( — ) particles in a 
NVT simulation at T = 0.7 and p = 0.004044 (in units of tr -3 ) as well as the aggregated number of McDonald's 
demon interventions per unit volume in a GCMD simulation with T = 0.7, S = 2.8658, and = 51 (• • •) 
over simulation time; bottom: pressure over simulation time for the NVT simulation ( — ) and the GCMD 
simulation (— ). 

Figure [2] Nucleus number per unit volume p n over nucleus size v from NVT simulation at T = 0.7 and p = 0.004044, 
with sampling intervals of 320 < t < 480 (o) and 970 < t < 1130 (o) after simulation onset, and from GCMD 
simulation with T = 0.7, S = 2.8658, and = 51 (•) in comparison with a prediction for the same conditions 
based on CNT (-). 

Figure \S\ Nucleation rate logarithm \a.J over supersaturation S at T = 0.65, 0.7, and 0.85 according to CNT (— ), 
the SPC modification ( — ) as well as HSL (■ • ■) compared to present GCMD simulation results (o). 

Figure [4] Nucleus temperature over nucleus size from GCMD simulation at T = 0.7 and S = 2.4958 for an intervention 
threshold size of = 15 (A), 30 (o), 48 (•), 65 (V), and 74 particles (■); dotted line: saturation temperature T a 
= 0.7965 of the vapor at constant pressure p = 0.134, which corresponds to the chosen supersaturation; dashed 
lines: guide to the eye. 

Figure [5] Intervention rate logarithm lnj7e over intervention threshold size of McDonald's demon during GCMD 
simulation at T — 0.7 and S — 2.4958 (□) in comparison with predictions based on CNT (— ) and the SPC 
modification ( — ); dotted line: CNT prediction shifted to the actual value of the nucleation rate; vertical line: 
critical size according to the SPC modification. 
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